Comparison script for SoilR model structures

1
jbeem@bgc-jena.mpg.de

Notes:

1. bulk & respired 14C used for contraints

• three time points for bulk 14C (2001, 2009, 2019)
• two time points for respired CO2 (2001, 2019)

2. steady-state C stocks (mean of 2001, 2009 stock data) used to determine inputs

This script requires the following input: obs.bulk.14c, obs.resp.14c, obs.Cstock

Workflow

1. Fit models by eye for initial parameter estimates
2. Run modfit using bulk and respired 14C time series data
3. Run bayesian parameter optimization
4. Reject models w/ stocks outside ±2SE of measured values

• Then what?

5. Adjust inputs w/ modeled stocks and reject models w/ unrealistic inputs
   • If inputs are not realistic [how to come up with realistic ranges?], may have to implement some constraints on parameter optimization, e.g. if (In > max | In < min) { error } else { fit pars}. Could calculate rejection rate associated with this process…
6. Compare models using AIC; other metrics? E.g. collinearity of parameters, ssr, bias, etc.
7. Get ages and transit times from “best fit” models, e.g. lowest AIC value
8. Compare distribution of ages/transit times with age distribution from density fractionation, thermal fractionation
9. Compare ages, transit times, to mineral characteristics (PCM, minC proportion, etc.)
10. Assess model structure/complexity [objective metrics here?]

Questions

1. How to set acceptable ranges for inputs?
2. How to deal with lag times? Would be nice to avoid fitting a lag parameter. Options: could set lag by depth, but likely the optimal lag value will vary by site.

Should copy relevant code for preparing obs.bulk.14c, obs.resp.14c

load pre-built bulk, respired 14C objs

SOC stocks

load min data

Error in names(test)[1, 3:9] : incorrect number of dimensions

Load initial parameter lists (1p, 2pp, 2ps, 3pp, 3ps)

modFit parameter optimization

Currently 3-pool lagged model fits are commented out. Fitting 6 parameters seems beyond the scope of the dataset. When the lag is not fit the modFit optimization performs better. Perhas could use the lags fit with the 2-pool model optimization for the 3-pool systems?

pool  3 < 0
pool  3 < 0
pool  3 < 0
pool  3 < 0
pool  3 < 0
pool  3 < 0
pool  3 < 0
pool  3 < 0
pool  3 < 0
pool  3 < 0

NaNs producedCannot estimate covariance; system is singularCannot estimate covariance; system is singularCannot estimate covariance; system is singularCannot estimate covariance; system is singularCannot estimate covariance; system is singularCannot estimate covariance; system is singularCannot estimate covariance; system is singularCannot estimate covariance; system is singularCannot estimate covariance; system is singularCannot estimate covariance; system is singularCannot estimate covariance; system is singularCannot estimate covariance; system is singularCannot estimate covariance; system is singularCannot estimate covariance; system is singularCannot estimate covariance; system is singularCannot estimate covariance; system is singularCannot estimate covariance; system is singularCannot estimate covariance; system is singularCannot estimate covariance; system is singularCannot estimate covariance; system is singularCannot estimate covariance; system is singularCannot estimate covariance; system is singularCannot estimate covariance; system is singularCannot estimate covariance; system is singularCannot estimate covariance; system is singularCannot estimate covariance; system is singularCannot estimate covariance; system is singularCannot estimate covariance; system is singularCannot estimate covariance; system is singularCannot estimate covariance; system is singularCannot estimate covariance; system is singularCannot estimate covariance; system is singularCannot estimate covariance; system is singularCannot estimate covariance; system is singularCannot estimate covariance; system is singularCannot estimate covariance; system is singularCannot estimate covariance; system is singularCannot estimate covariance; system is singularCannot estimate covariance; system is singularCannot estimate covariance; system is singularCannot estimate covariance; system is singularCannot estimate covariance; system is singularCannot estimate covariance; system is singularCannot estimate covariance; system is singularCannot estimate covariance; system is singularCannot estimate covariance; system is singularCannot estimate covariance; system is singularCannot estimate covariance; system is singularCannot estimate covariance; system is singularCannot estimate covariance; system is singularCannot estimate covariance; system is singularCannot estimate covariance; system is singularCannot estimate covariance; system is singularCannot estimate covariance; system is singularCannot estimate covariance; system is singularCannot estimate covariance; system is singularCannot estimate covariance; system is singularCannot estimate covariance; system is singularCannot estimate covariance; system is singularCannot estimate covariance; system is singular

modFit timeseries

modFit SSR

modFit AIC

Lagged 2-pool models have higher AIC values than the non-lagged models for all sites except for the deeper layers at the GRrf sites (for which the models fit poorly anyway).

When comparing lagged 1-pool models to the 2-pool models without a lag, the 2-pool models perform better at the soil surface, but the 1-pool models perform slightly better at the deepest depth. Performance for the 10-20 cm layer is mixed. In general, the 1-pool models perform better than 2-pool models for the warm climate sites, while the 2-pool models perform better at the cold climate sites. For the cool climate sites, the 2-pool models perform better for andesitic parent materials; however, basaltic and granitic soils follow the overall depth trends: with better performance of the 2-pool models at the soil surface, and better performance of the 1-pool models at the deepest depth.

modFit par sens

Sensitivity of parameter kslow is 0! 

modFit identifiability

Parameter combination of length > 2 are rarely identifiable

Inputs

For the 2pp models the sites ANpp 0-10 and BSwf 10-20 & 20-30 had much higher inputs than the other sites. For the 2ps models it’s the same, with the addition of GRwf 10-20. Sites ANpp 0-10 and BSwf 20-30 also have anomalously high mean ages. All of the sites with high inputs have transit times ~ 1 yr (TT for ANpp 0-10 = 9 yr), in contrast to values ranging from 60 yrs to hundreds of yrs.

Inputs for 3pp models not much different from 2pp models, alhough there is a tendency to be slightly greater.

Call:
lm(formula = age ~ pm, data = .)

Residuals:
    Min      1Q  Median      3Q     Max 
-172.93  -67.50   20.40   38.09  182.54 

Coefficients:
            Estimate Std. Error t value  Pr(>|t|)    
(Intercept)   329.92      57.56   5.732 0.0000692 ***
pmbasalt     -168.55      74.31  -2.268    0.0410 *  
pmgranite    -179.56      74.31  -2.416    0.0311 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 115.1 on 13 degrees of freedom
Multiple R-squared:  0.3462,    Adjusted R-squared:  0.2457 
F-statistic: 3.443 on 2 and 13 DF,  p-value: 0.06312


Call:
lm(formula = age ~ pm * eco + middepth, data = .)

Residuals:
    Min      1Q  Median      3Q     Max 
-262.00  -46.52   -9.49   24.21  814.72 

Coefficients:
                  Estimate Std. Error t value   Pr(>|t|)    
(Intercept)        499.635     88.795   5.627 0.00000185 ***
pmbasalt          -381.732    103.731  -3.680    0.00072 ***
pmgranite         -490.050    103.731  -4.724 0.00003130 ***
ecocool            241.939    103.731   2.332    0.02508 *  
ecocold           -202.389    103.731  -1.951    0.05845 .  
middepth             6.187      3.336   1.854    0.07145 .  
pmbasalt:ecocool  -199.109    156.489  -1.272    0.21098    
pmgranite:ecocool  -42.803    146.698  -0.292    0.77205    
pmbasalt:ecocold   240.108    146.698   1.637    0.10994    
pmgranite:ecocold  340.756    182.739   1.865    0.06996 .  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 179.7 on 38 degrees of freedom
Multiple R-squared:  0.676, Adjusted R-squared:  0.5992 
F-statistic: 8.808 on 9 and 38 DF,  p-value: 0.0000005264

Joining, by = "PMeco"

Call:
lm(formula = tt ~ pm * MAT + middepth, data = .)

Residuals:
     Min       1Q   Median       3Q      Max 
-109.975  -24.356    6.488   29.804  184.610 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)    172.767     54.457   3.173 0.002862 ** 
pmbasalt      -281.069     69.666  -4.035 0.000233 ***
pmgranite      201.070    119.927   1.677 0.101228    
MAT            -11.119      5.842  -1.903 0.064061 .  
middepth         9.410      1.009   9.323  1.1e-11 ***
pmbasalt:MAT    20.567      7.442   2.764 0.008518 ** 
pmgranite:MAT  -27.598     12.463  -2.214 0.032416 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 55.73 on 41 degrees of freedom
Multiple R-squared:  0.7688,    Adjusted R-squared:  0.735 
F-statistic: 22.72 on 6 and 41 DF,  p-value: 1.344e-11

Ages & transit times

Fits for the 2pp and 2ps models were similar across the 0-10 cm depth layer. The only noticable exceptions were for the GRwf and GRrf sites. For the GRrf site, the fit for the 2pp model was essentially a one pool model, with nearly identical tau values of 1/170 for the fast and slow pools. In contrast, tau for the fast pool in the 2ps model fit at the GRrf site was 1, while the slow pool tau was 1/170. For the GRwf site, the 2pp fit had tau values of 21 and 152 for fast and slow pools, respectively, while the 2ps taus were smaller (longer turnover time): 57 and 1075 respectively.

In general, due to the similarities in the fits and the differences in model structure, the ages and transit times were substantially longer for the series models than for the parallel models. However, mean ages and transit times for both the 2pp and 2ps fits appeared to be excessively long at several of the sites. This problem was especially apparent for the 10-20 cm and 20-30 cm depth layers, although we also observed ages for the ANpp site in the 0-10 cm depth layer that exceeded those of the other sites by an order of magnitude (ca. 1500 vs. 150). For the 10-20 cm layer, the cold climate granite site had the oldest mean ages: 1750 and 2500 for the 2pp and 2ps models respectively, with the mean age of the remaining sites = 350 yrs. The most striking outlier for age was the 20-30 cm layer at the BSwf site, for which both the 2pp and 2ps models had mean ages ≥ 95,000, exceeding the mean ages of the other sites by 1 to 2 orders of magnitude.

Bayesian fits

MCMC chains for most of the parameters, across all of the model structures, fail to converge after 1000 iterations. Should try 5000 itr for a subset, evaluate improvement, and consider optimal number of iterations.

Can also consider using smaller updatecov value (currently set equal to the number of iterations) during the burnin period to harness the power of the adaptive Metropolis method. Delayed rejection algorithm is probably not necessary given the generally high level of accepted runs (most >90%).

MCMC AIC

Comparison of AIC values between modFit parameter sets and MCMC sets indicate a clear performance increase for the MCMC parameters. However, when comparing relative AIC values among the MCMC parameter sets fit for parallel and series model structures, and the structures with or without lags are essentially indistinguishable. The difference in AIC values, i.e. \(\Delta\)AIC, is within 2 for all of the sites, indicating a high level of support for either model structure and either parameter set. Input calculations are also similar across model structures, although more of the two pool series model fits yield unrealistically high (or low) inputs.

[1] 0.5945797

[[1]]

[[2]]

[[3]]

[[4]]

[[5]]

[[6]]

[[7]]

[[8]]

[[9]]

[[10]]

[[11]]

[[12]]

[[13]]

[[14]]

[[15]]

[[16]]

[[17]]

[[18]]

[[19]]

[[20]]

[[21]]

[[22]]

[[23]]

[[24]]

[[25]]

[[26]]

[[27]]

MCMC Inputs

Results from the inputs suggest that some of the parameter sets are simply not realistic. The question remains as to whether or not the remaining inputs are realistic.

There is no discernable trend with climate, but there does appear to be a trend with depth. Surprisingly, perhaps, the 0-10 cm depth does not diverge substantially from the 10-20 cm depth (I would expect the surface inputs to be substantially higher, primarily due to root density).

Furthermore, the andesitic soils appear to have higher inputs than the other soils. This could be true, but I don’t know if I believe it. I wonder if there is a way to demonstrate parent material effects on GPP, e.g. as a function of nutrient status? Likely GPP-based estimates would be sufficient, as I wouldn’t expect partitioning between above and belowground inputs to be affected by parent material. Craig has mentioned that if there were to be any impacts of parent material on fertility, and by extension on belowground inputs, the basalt soils would likely be the most fertile due to inherited phosphorus.

Error in lapply(SA.TT.2pp.lag.ls, "[[", "ins") : 
  object 'SA.TT.2pp.lag.ls' not found

Toss list: ANpp 0-10: 2ps, lag BSwf 10-20: 2ps, lag; 2pp, lag BSwf 20-30: 2ps, lag; 2pp, lag GRwf 10-20: 2ps, lag

MCMC SOC stocks

MCMC SA & TT

Ages

System ages for the series model fits tend to be older than the parallel fits. Additionally, some of the ages are unrealistically old: ANpp 0-10 (32744). Ages for both the parallel and series fits are extraordinarily old for the BSwf_20-30 site, although the series model age is an order of magnitude greater than the parallel model: series = 921536 y, parallel = 94760 y. The age for the ANpp_0-10 site is also too old for parallel model: 1032 y (although still an order of magnitude younger than the series model age for the same site and depth).

Transit times

Transit times for the parallel models seem mostly reasonable. The cool climate basalt soils from the two deepest layers are an exception in that they are nearly 0, indicating a failure on the part of the model or poor data. Specifically, the reason that BSwf_20-30 is the worst fitting site for both model structures is likely related to a failure to resample the same site, due to logging in 2019 and possibly a failure to locate the correct site in 2009.

In contrast, the two deepest GRrf sites have overly long transit times for both model structures, likely a function of the anomolously old respired 14C signal at this site for these depths. The trend is also potentially apparent in the surface soil, where the transit time is nearly 150 years for the parallel model, compared to 75 years for the andesitic and basaltic soils. However, this is not visible in the series model for this site and depth.

Overall, the transit times for the series models do not seem realistic. Perhaps this is reflected in the estimates for inputs as well? Would be interesting to see if the lagged models perform better.

SA & TT for lagged models

The lagged models show similar trends to the models without lags, with more sites returning both ages and transit times that are unrealistically old for the series model structure than the parallel structure.

Oddly I don’t observe large differences in either age or transit time between the lagged and unlagged model fits. The exceptions are the sites where the model appears to have found local minima, and the age/transit time statistic is unreasonably high (or low) for either the lagged or unlagged fit.